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Abstract 

Titan's moment of inertia, calculated assuming hydrostatic equilibrium from gravity field data obtained during 
the Cassini-Huygens mission, implies an internal mass distribution that may be incompatible with complete differ- 
entiation. This suggests that Titan may have a mixed ice/rock core, possibly consistent with slow accretion in a 
gas-starved disk, which may initially spare Titan from widespread ice melting and subsequent differentiation. A 
partially differentiated Titan, however, must still efficiently remove radiogenic heat over geologic time. We argue 
that compositional heterogeneity in the major Saturnian satellites indicates that Titan formed from planetesimals 
with disparate densities. The resulting compositional anomalies would quickly redistribute to form a vertical density 
gradient that would oppose thermal convection. We use elements of the theory of double- diffusive convection to 
create a parameterized model for the thermal evolution of ice/rock mixtures with a stabilizing compositional gradi- 
ent. Simulations are performed for a wide range of initial conditions to account for large uncertainties in material 
properties and accretionary processes. Ultimately, for realistic density gradients, double-diffusive convection in the 
ice/rock interior can delay, but not prevent, ice melting and differentiation, even if a substantial fraction of potassium 
is leached from the rock component. Consequently, Titan is not partially differentiated. 
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1. Introduction 

Titan is Saturn's largest satellite, the second largest 
in our Solar System. The surface of Titan features sta- 
ble liquid methane/ethane lakes, a hydrocarbon precip- 
itation cycle comparable to Earth's hydrology, and myr- 
iad additional geologic features such as dune fields and 
mountainous terrain (Jaumann et al. , 2009). From an 



astrobiology perspective, Titan is interesting as a mod- 
ern analogue for pre-biotic Earth (|Raulin et al. , 2009). 



Titan has a thick nitrogen atmosphere, which must 
be sustained by the continuous production of methane 
(Atreya et al. , 2006). Methane may be stored subsur- 



face or perhaps more deeply in stable clathrate-hydrates 
along with other volatiles and episodically outgassed 



when internal heating causes melting ( Gautier and Her- 



sant 



2005 Tobie et al. , 2006). Understanding Titan's 



fascinating surface, atmosphere, and chemistry requires 
knowledge of the evolution of its interior. 

Measurements of the gravity field of Titan can place 
indirect constraints on the current structure of its in- 
terior. Titan's gravity harmonics were determined to 
degree 3 through careful tracking of the Cassini space- 
craft during four flybys (less et al. 2010). The calcu- 



lated ratio J2/C22 is consistent with the value of 10/3 
that is itself consistent with a gravity field dominated 
by a nearly hydrostatic quadrupole, although non-zero 
values of other degree 2 and 3 coefficients indicate that 
nonhydrostatic features are present. If hydrostatic equi- 



librium is assumed, however, then Titan's moment of 
inertia (Mol) factor is found to be C ~ 0.34. Here, the 
Mol is CM S R 2 S1 where M s and R s respectively represent 
the mass and radius of Titan. Measurement of the tidal 
Love number k<i from additional flybys reveals a rela- 
tively large response of the gravity field to the Saturnian 
tidal field, indicating the presence of a global, subsur- 
face ocean ( |Iess et aL 2012). The decoupling of Titan's 
shell from its interior with an ocean may also explain 



Titan's long- wavelength topography (Nimmo and Bills 
2010), spin pole orientation (Bills and Nimmo, 2011), 



and Schumann resonance ( |Simoes et al. , 2012). 



The 2-<j error on J2/C22 is about 3%, centered on 



10/3 and dominated by the error in J2, according to less 
et al. (2010). We can write J2 = J2,h+J2,nh and likewise 



C22 = C22,h J r C 22, nh-> where "h" means hydrostatic and 
"nh" means non- hydrostatic, so J2,h/C22,h is exactly 
10/3. If the non- hydrostatic parts were completely un- 
correlated, then the closeness of the observed value to 
10/3 is highly significant and precludes a substantial er- 
ror in Mol. If, on the other hand, the non-hydrostatic 
parts tend to be correlated (as they would be, for ex- 
ample, in the tidal heating model of |Nimmo and Bills 
(2010) or if Titan has undergone True Polar Wander to 



the preferred orientation of the non-hydrostatic part of 
the Mol), then the small deviation away from 10/3 does 
not guarantee smallness of the non-hydrostatic part and 
the Mol is accordingly uncertain. Given the observed 
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power in degree 3 gravity, less et al.| ( |2Q1Q[ ) propose that 
the Mol could be as low as 0.33 and we adopt this as 
a plausible lower bound. The true Mol is likely to be 
smaller than 0.34, the value inferred from J2/C22, be- 
cause this is a lower rotational energy state. Since a 
fully differentiated Titan may allow 0.33 (but perhaps 
not 0.34), the differentiation state of Titan is accord- 
ingly not yet determined from observation. 

Titan's Mol can be compared to previous results re- 
garding other icy satellites. In particular, Titan's Mol 
coefficient is intermediate to the previously-measured C 
~ 0.31 for Ganymede (Anderson et al. 1996a) and C ~ 
0.35 for Callisto (Anderson et al. 2001), where hydro- 



static equilibrium was assumed a priori for the Galilean 
satellites. While Ganymede can easily be modeled as 
a differentiated satellite with an iron core under a rock 



shell and an outer ice layer (e.g., Sohl et al. , 20021), mod- 
els of a differentiated Callisto are not consistent with the 



reported Mol (e.g., Anderson et al. , 2001). A proposed 
interior structure for Callisto features a rock/ice core 
with rock mass fraction near the close packing limit, 
with an overlying icy mantle that was depleted of rock 
by Stokes settling (Nagel et al. , 2004). In this case, 
an ice/rock lithosphere in which density decreases with 
depth overlays the icy mantle. Titan might represent an 
intermediate case between differentiated Ganymede and 
undifferentiated Callisto. However, the inferred (and 
widely cited) partially differentiated state for Callisto 
is based on a gravity inversion that assumes J2/C22 = 
10/3, and perhaps the inferred Mol is incorrect. As 
for Titan, the sense of the error is likely to cause an 
overestimate of the Mol because non-hydrostatic con- 
tributions to J 2 and C22 are likely to be positive. At 
present, we cannot exclude the possibility that all large 
icy satellites are fully differentiated. (In this paper, we 
use that term to imply the complete separation of ice 
from rock; the further differentiation of an iron core 
from the rock is a separate issue that we do not address.) 
It should be noted that non-hydrostatic effects are more 
likely to be important in slowly rotating bodies (Titan 
and Callisto) relative to Ganymede because the hydro- 
static effects scale as rotational frequency squared for 



both tides and synchronous rotation (Gao and Steven- 



son 



2012, and in preparation). 



Certainly Callisto and Ganymede are different in ap- 
pearance and Callisto, like Titan, lacks a global mag- 
netic field. Relative differences in the experienced in- 
tensity of the Late Heavy Bombardment (LHB) may 
explain the Ganymede/Callisto dichotomy. In the so- 
called "Nice model," the gas giant planets swiftly re- 
aligned once Jupiter and Saturn crossed their mean 
motion resonance ~700 million years after planetary 
formation, quickly causing the outward migration of 
Uranus and Neptune and the evolution of Jupiter and 



Saturn's orbital eccentricities (Tsiganis et al. , 2005). 
The ensuing destabilization of the planetesimal disk and 
the asteroid belt caused the LHB of both the outer and 



inner solar system rtGomes et al. , 2005). Ganymede is 



siderably more high-energy impacts during the LHB. 
For assumptions about the planetesimal disk consistent 
with the Nice model, the differences in received energy 
during the LHB are sufficient to cause Ganymede to 
differentiate, while Callisto's rock and ice components 



may remain unmixed (Barr and Canup, 2010). If Titan 



survived both accretion and the LHB without melting 
the ice in its deep interior, then it could have remained 
partially differentiated like Callisto, at least initially. 

Gas giant satellites like Titan were accreted from the 
outskirts of the disk of material surrounding their par- 
ent planets. Nascent gas giant planets must accrete 
enough rock such that they can also accumulate large 
amounts of gas, principally hydrogen and helium, before 
the dissipation of the protoplanetary disk, which had a 
lifetime of a few million years (e.g., Lissauer and Stevem] 



son 



2007 Lunine et al. , 2009 ) . In the core- nucleated gas 



accretion model, for instance, the formation of Saturn 
began with the accretion of kilometer-sized rocky plan- 
etesimals from the minimum mass sub nebula (MMSN) 



(e.g., Lissauer and Stevenson 2007). However, if gas 



giant satellites were accreted from the outskirts of a 
disk of planet-forming material as dense and gas-rich 
as the MMSN, then Jupiter's Galilean satellites (and, 
by analogy, Titan) would have formed quickly and hot 
(e.g., |Mosqueira and Estrada| |2003a| ) and therefore dif- 
ferentiated. 

Internal structure models for Titan often assume 
widespread ice melting and thus differentiation. Many 
studies investigated the thermal evolution of Titan as- 



suming a large silicate/iron core (e.g., Sohl et al. , 2003 



Tobie et al.| [2005]), before Mol data from |Iess et al. 



(2010) cast doubt on such models. Convection in a sil- 



icate core has been invoked to melt clathrate hydrates 



and cause episodic outgassing of methane (Tobie et al 



2006 ) . Another class of models assumes that Titan con- 



tains a large core of hydrated silicates, chiefly the ser- 



pentine mineral antigorite ( Grindrod et al. 



2008 



eTaTl [20071 ICastillo-Rogez and Lunine| |2010 



2012 



Fortes 



ortes. 



Tobie et al.[ 2012). Serpentinization of silicates 



in icy satellites is likely to be rapid in the presence of 
liquid water, e.g. during differentiation or subsequent 



hydrothermal convection (Ransford et al. , 1981 Travis 
|and Schubert 2005 ) , because the serpentinization reac- 



tion promotes crack formation and increasing material 
permeability (MacDonald and Fyfe, 1985[ ). An inte- 
rior dominated by hydrated silicates is consistent with 



constraints on Titan's chemical evolution (e.g., Fortes 



et al., 2007 Fortes, 2012). But a partially differenti- 



ated Titan, in which Titan's deep interior is a mixture 
of ice and rock, is still consistent with the gravity data 



closer to Jupiter than Callisto, so it likely suffered con- 



( |Iess et al. , 2010), and such models have not been fully 
vetted. 

The accretion of undifferentiated icy satellites is often 
assumed to have occurred in a "gas-starved" disk, where 
material similar to that found in the gas giant's feed- 
ing zone is fed to the accretion process over millions of 
years. The resulting long satellite accretion timescales 
allow many of the Galilean satellites, particularly Cal- 



2 



listo, to avoid complete differentiation (e.g., Canup and 



Ward, 2002). In another model, the major satellites 
of Jupiter and Saturn formed in a solid-enhanced min- 
imum mass planetary nebula, avoiding differentiation 
as satellites open gaps in the nebula ( |Mosqueira and| 
Estrada, 2003a b). With a particular compositional gra- 
dient in the initial circumplanetary disc, the composi- 
tions of the major Saturnian satellites can be roughly 



reproduced (Mosqueira et al. , 2010). In any model, 



accretion must also be delayed to escape intense, but 
short-lived, radiogenic heating from the decay of 26 Al 
if an icy satellite is to avoid complete differentiation 



(e.g., |McKinnon[ |1997| |Barr and Canup| |2008| ). De- 
spite myriad threats to unmelted ice, slow accretion in 
a gas-starved disk permits the formation of a partially 
differentiated Titan (iBarr et all 120101). 



The purpose of this study is to determine whether 
a partially differentiated Titan is stable over geologic 
time. After accretion, the internal ice/rock mixture 
must efficiently expel radiogenic heat, which would oth- 
erwise melt the ice and allow irreversible sinking and 
separation of the rock component. Because the major 
Saturnian satellites are remarkably heterogenous, we ar- 
gue that Titan likely accreted from planetesimals with 
disparate rock mass fractions. Even without that het- 
erogeneity, there may be a tendency to form a stably 
stratified interior at the outset. Soon after accretion, 
a stabilizing density gradient would be established in 
which rock mass fraction increases with depth. As- 
pects of the theory of double-diffusive convection allow 
us to formulate a one-dimensional, parameterized ther- 
mal evolution model in which a convecting layer slowly 
grows to encompass the entire ice/rock interior. We 
perform simulations for a large set of initial conditions 
because of uncertainties in key parameters and to an- 
alyze the sensitivity of our results. Finally, we discuss 
alternatives to a partially differentiated Titan and con- 
sider possible implications for other icy satellites. 



2. Theoretical formulation 

Our simple model assumes that Titan is a sphere of 
radius R s and mass M s , consisting of an undifferenti- 
ated ice/rock interior with radius Ri and an overlying 
shell. The structure of the outer shell may be compli- 
cated, especially considering the likely presence of an 
ocean, but we do not need to model it in detail. We 
perform parametrized thermal evolution simulations to 
assess whether Titan can remain partially differenti- 
ated until the present. Stabilizing compositional gradi- 
ents, caused by increasing rock mass fraction and thus 
density with depth, likely inhibit convection in Titan's 
ice/rock interior. Elements from the theory of double- 
diffusive convection are incorporated to model thermal 
evolution with this complication. Although this simple 
model will not precisely reconstruct the thermal his- 
tory of Titan, it does provide a first-order test of the 
hypothesis of partial differentiation. 



2.1. Establishment of vertical density gradients 

The major Saturnian satellites exhibit significant 
compositional heterogeneity. The mass-weighted den- 
sity of the major Saturnian satellites, excluding Ti- 
tan (Mimas, Enceladus, Tethys, Dione, Rhea, Hype- 
rion, Iapetus, and Phoebe), is p av ~ 1,270 kg m -3 
flJacobson||2004l ). The average densities of each satel- 



lite vary widely, from less than 1,000 kg m -3 to over 
1,300 kg m -3 . In fact, the inner satellites (the first five 
of the aforementioned) and Saturn's almost-pure water 
ice rings may have been formed from tidal stripping of a 
Titan-sized satellite that was lost within Saturn's clas- 



sical Roche limit (Canup, 2010). Titan's mean density 
is ~1,881 kg m -3 , but its uncompressed density is not 
much higher than the densities of the other Saturnian 
satellites. 

Since Titan was plausibly formed from planetesimals 
of size and composition similar to these satellites, lat- 
eral density anomalies as large as 10% likely existed 
after accretion. These unstable anomalies drive a flow 
that will restructure the ice/rock mixtures so that more 
dense mixtures underlie less dense mixtures, in a length 
of time that scales as 



tf 



ApiatgD 



(i) 



where /i is viscosity, Api at is the horizontal density 
anomaly, g is gravitational acceleration, and D is the 
length scale of the anomalies. A competing process is 
the Stokes settling of the dense rock fragments from the 
ice/rock mixture. Using the usual Stokes settling veloc- 
ity, the time in which the rock fragments settle out of 
the mixture may be written as 



ApgR 2 ' 



(2) 



where Ap is the density difference between the ice and 
the rock (on the order of 1,000 to 2,000 kg m -3 ) and 
R is the radius of the rock fragments. For R ~ 1 m 
and D ~ 100 km, the formation of the vertical density 
gradients overwhelms the Stokes settling of the rocks 
if the lateral density anomaly is greater than only one 
part in ten billion of the density difference between the 
ice and rock, i.e. tf < t s if Api at > 10 _10 Ap. Clearly, 
this is easily achieved. 

Mixing during accretion will not homogenize the dif- 
ferences arising from large planetesimals. Importantly, 
lateral differences during the accretion lead to vertical 
differences in ice/rock ratio after settling, even if indi- 
vidual impact events mix vertically. Therefore, models 
of Titan's thermal evolution must consider the existence 
of up to ~10% stabilizing density gradients in Titan's 
ice/rock interior initially (and even smaller effects are 
still very important because thermal expansion of ice 
produces much smaller density differences). A vertical 
density gradient, where density increases with depth, 
opposes thermal buoyancy, inhibiting convection and 
thus preventing substantial heat transfer. If radiogenic 
heat is not removed from the system, then the ice in 



3 



Titan's deep interior will melt, especially if low melting 
point components such as ammonia are present. The 
liberated rock components will sink as the melt per- 
colates upwards. Runaway differentiation can occur if 
the gravitational energy released by the sinking of the 



rock component is sufficient to melt more ice (Fried 



son and Stevenson 1983), but complete differentiation 



may occur without invoking runaway differentiation if 
convection is sufficiently suppressed. 

2.2. Double- diffusive convection applied to Titan 

The theory of double-diffusive convection was orig- 
inally developed to elucidate various phenomena in 
oceanography (Stern, 1960). In general, double- 



diffusive phenomena occur in the presence of two op- 
posing density gradients, where the two responsible 
components have different molecular diffusivities (e.g., 



Turner, 1974). The theory has been extended to the dy- 



namics of magma chambers and stellar interiors (e.g., 



Spiegel) |1972| |Schmitt[ |1983| |Turner| |1985| ), but it re- 
mains most often utilized in oceanography. For in- 
stance, in Earth's high latitudes, the polar oceans are 
composed of cold, fresh water overlying warm, salty wa- 
ter (e.g., |Schmitt| |1994[ ). In the presence of the stabi- 
lizing salinity gradient, convective motions resembling 
Rayleigh-Benard convection arise after a stage of steady 



oscillating motions in layers cooled from above (Vero- 
nis , 1968} jNoguchi and Niino} |2Q1Q| ) . These well-mixed, 



convecting layers increase in thickness and frequently 
merge with neighboring layers of equal density. This is 
known as the "diffusive-layer" mode, in which salinity, 
the stabilizing component, has a much lower molecu- 
lar diffusivity than temperature, the component driving 
convection. 

Titan's deep interior features analogous opposing 
temperature and compositional gradients. The compo- 
sitional gradient is simply a measure of the increasing 
rock mass fraction in the ice/rock mixture with depth. 
Although a temperature gradient can induce particle 



motion through interfacial premelting (Rempel et al. 



2001), the relevant time scales are far too long for ap- 



plication to icy satellites. Hence, the molecular diffu- 
sivity associated with the rock component is effectively 
zero. Because the molecular diffusivity of salt is non- 
zero, in contrast, complicated temperature and salinity 



profiles exist in the ocean (e.g., |Huppert and Turner 



1981), with time-dependent fluxes of salt and temper- 



ature across diffusive interfaces between the multiple 
convecting layers (e.g., |Worster| [2004). Regardless, the 
term "double-diffusive" is still applicable to our model 
of convection in the presence of a stabilizing gradient. 
As outlined below, aspects of the theory applied to Ti- 
tan do not depend on the molecular diffusivity of the 
component opposing temperature. 

In Titan, double-diffusive convection occurs in the 
presence of a stabilizing rock mass fraction gradient as 
its deep interior is cooled from above or heated from 
within. The density within the ice/rock interior is a 



linear function of rock mass fraction and temperature 



flTurner||1974| ): 



p = p (l-aAT + f3AS), 



(3) 



where po is the density at some reference depth, a is the 
coefficient of thermal expansion, and f3 = (1/ p)(dp/dS) 
is the coefficient of compositional expansion, analogous 
to a, associated with the change in rock mass fraction. 
The change of density due to pressure alone can be omit- 
ted because it does not contribute to the convective sta- 
bility and is in any event too small to greatly affect the 
form of the convective motions. Table [T] contains nu- 
merical values of these coefficients and other constants 
used in this study. The rock mass fraction may be cal- 
culated ( |Barr et al.[|2010| ): 



Pr(p~ Pi) 



(4) 



P(pr-PiY 

where p is the average density, p r is the rock density, and 
pi is the ice density. High-pressure ice phases V, VI, and 
VII exist within Titan's deep interior; a represent ative 
density pi = 1,400 kg m -3 is chosen in accord with|Barr] 



et al. (2010). Although the mineralogy of Titan's rock 



component is not completely known, CI carbonaceous 
chondrites and Prinn-Fegley rock, an assemblage that 
models condensate from the proto- Jovian nebula, are 
possi ble analogues, motivating the choice p r = 3,300 k g 
m~ 3 ( |Mueller and McKinnon} 1 1988} |Barr et al.[|2010 ). 

A quick calculation with Eqs. [3] and [4] reveals why 
compositional gradients make melting much more likely 
in a partially differentiated Titan. A ~10% density con- 
trast across Titan's interior corresponds to AS ~ 10 _1 
between ice/rock parcels at the top and bottom of the 
mixture. So, an (impossible) temperature increase of 
AT rsj 10 3 K is required before compositionally dense 
parcels can become buoyant. 

For cooling from above or for steady heating from 
within, but with the top temperature kept constant, 
convection begins at the top of Titan's ice/rock interior. 
With a linear density distribution, an equation for a 
measure of the initial stability gradient may be written 



(Turner 1973): 



dz ' 



(5) 



where Ns is a frequency of oscillation. When heat flows 
out of the top of the convecting layer, the associated 
buoyancy flux may be calculated (Turner, 1973): 

gaH 



B = 



pC p 



(6) 



As cooling continues from above, a well-mixed convect- 



ing layer grows with height (rTurner, 1973): 



h 



2Bt\ 



1/2 



(7) 



Therefore, the rate of the growth of the convecting layer 
may be calculated: 

1/2 

(8) 



dh 



aH 



_2pC p P(dS/dz)t 
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Constant 


Definition 


Value 


Units 


Ref. 


A 


Activation parameter for ice 


25 




[1] 


K 


Thermal diffusivity 


10" 6 


9 —1 

m s 


[1] 


a 


Coemcient ol thermal expansivity 


10 4 


-rr — 1 

K 1 


r-i i 

[i] 


n 
P 


Coefficient of compositional expansivity 


1.1 


- 




M s 


Mass of Titan 


1.345 x 10^ 


kg 


[2] 


T) 

R S 


Kadms oi litan 


ZOYO 


km 




Ri 


Radius of ice/rock interior 


2000 


km 


* 


Pi 


Density of ice 


1400 


kg m -3 


[2] 


Pr 


Density of rock 


3300 


kg m -3 


[2] 




Specific heat of ice 


2100 


J kg" 1 K- 1 


[2] 




Specific heat of rock 


900 


J kg" 1 K- 1 


[3] 



Table 1: List of key model constants. References: 1. Friedson and Stevenson (1983), 2. Barr et al. (2010j, 3. |Grindrod et al.| ( |2008| >, 
*this study. 



Water/Ice Shell 



H' 



B. L 

Ice/Rock 
Convection 

Ice/Rock 
Interior 



where H is the heat flux out of the convecting ice/rock 
layer into the overlying shell, Q is radiogenic heat pro- 
duction, fi is the volume of melted ice, and L m is the 
latent heat of melting. The specific heat of the ice/rock 
mixture is a mass-weighted average ( |Barr and Canup| 
[2008] ): 

(10) 



Cp — SCp^ r + (1 — S)Cp^i, 



where C Pyr and C p j are the specific heats of the rock 
and ice components, respectively. Assuming that no 
melting occurs, we can conveniently rewrite Eq. [9j 



m 

dt 



3H iff 



(ii) 



Figure 1: Cartoon showing the assumed structure of Titan af- 
ter the initialization of double-diffusive convection in the interior 
ice/rock region, before the convecting layer has grown to encom- 
pass the entire interior. Dimensions are not to scale and the 
overlying water/ice shell is not drawn in detail. Key model pa- 
rameters are indicated as used in the text. 



Figure [T] shows the internal structure of Titan af- 
ter double-diffusive convection has proceeded for some 
time. The radius of the ice/rock interior is Ri. The non- 
convecting ice/rock core has radius R C1 where T c and T m 
respectively represent the outer and inner temperatures. 
A thermal boundary layer exists between the convect- 
ing outer shell and the convecting ice/rock interior with 
thickness S. The temperature at the top of the thermal 
boundary layer, T&, is held constant, while the tempera- 
ture of the convecting layer is . Parametrized thermal 
evolution models typically assume adiabatic tempera- 
ture gradients in convecting layers, but these tempera- 
ture increases with depth are negligible. 



2.3. Governing equations 

The governing equation is the conservation of energy 
equation for the convecting ice/rock layer: 

4nR*H = ^p(Rl - Rl) (q - C p ^\ - p^L^ (9) 



We consider contributions to radiogenic heat produc- 
tion from 40 K, 235 U, 238 U, and 232 Th. Radiogenic heat- 



ing in the ice/rock mixture may be calculated (e.g., Barr 
and Canup||2008| ): 



Q(t) 



£^c n;0 P n; oexp(- 



(12) 



where, for the nth isotope, c n ^ is the initial abundance, 
P n $ is the initial specific heat production, and A n is 
the decay constant. Table [2] lists the constants used 
to calculate radiogenic heat production. Because con- 
tact between silicates and ammonia/water liquid can 
leach volatile 40 K out of the rock component (e.g 



En- 



gel et al. , 1994), we multiply the initial abundance of 



potassium by /, a depletion factor. Figure [2] shows rep- 
resentative temperature increases over 4.5 Gyr for S = 
1.0 and 0.8 from radiogenic heat production with / = 
1.0, 0.7, and 0.3. 

2.4- Thermal evolution of Titan's ice/rock interior 

The thermal evolution the ice/rock mixture within 
Titan may be simulated with a one-dimensional 
parametrized model. A critical parameter, subject to 
large uncertainty, is the viscosity of the ice/rock mix- 
ture. Although the deformation of ice is non-Newtonian 
under the conditions of Titan's deep interior, it can be 



5 



Isotope 


C n ,0 [PPb] 




[i/^y r J 


238jj 


16.91 


93.7 


0.155 


235|j 


5.275 


569 


0.985 


232^ 


38.71 


26.9 


0.0495 


40 K ^ 40 Ar 


813.8 


1.02 


0.0581 


40 K 40 Ca 


813.8 


26.69 


0.4962 



Table 2: Constants used to calculate radiogenic heat production, adapted from|Barr and Canup| j2008| 
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Figure 2: Temperature increase in a parcel of rock (black, S = 
1.0) and in an ice/rock mixture with S = 0.8 (red). Solid, dashed, 
and dotted lines respectively represent / = 1.0, 0.7, and 0.3. 



modeled with an Arrhenius-like equation for low strain 



rates (Friedson and Stevenson 1983): 



exp 



A 



T 



(13) 



where A is an activation parameter and T m is the melt- 
ing point of water ice, calculated using data from |Pe-| 
jtrenko and Whitworthj ([1999 ) . For conditions appropri- 
ate to icy satellites, adding rock particles to a water ice 
matrix may cause viscosity to increase by an order of 
magnitude (|Friedson and Stevenson, 1983), but we do 



not explicitly model these effects in this study. 

At first, a conductive boundary layer grows with 
thickness 5 ~ \f^t at the top of the ice/rock interior. 
The conductive heat flux out of the boundary layer is 



H 



kAT 



(14) 



where k = tzpC p is thermal conductivity and k is ther- 
mal diffusivity. The initial growth of the conductive 
boundary layer stops and convection begins when the 
Rayleigh number for the boundary layer reaches a crit- 
ical value, i.e., 



V(T)K 



(15) 



where AT is the temperature contrast across the bound- 
ary layer and S c is the critical thickness of the boundary 
layer. 



After convection begins, the viscosity contrast be- 
tween the convecting ice/rock mixture and the overlying 
shell of high pressure ice is less than the ~4 orders of 
magnitude required to produce a stagnant lid. So, the 



heat flux is given by the following equation ( jFriedson 



and Stevenson 1983): 



H = O.lk(AT) 



4/3 



got 



1/3 



(16) 



where AT = T5 — T^ is the temperature contrast driv- 
ing convection. The thickness of the thermal boundary 
layer, <$, can be calculated using the condition for con- 
vective instability 



gaATSf 



rj(Ti - 0.5AT)*; 



= 10 6 



(17) 



In reality, the thermal evolution of the ice shell above 
this thermal boundary layer is complicated. Titan likely 
has an shell consisting of an outer ice I layer, an ocean, 
and a mantle of high-pressure ice polymorphs (e.g.,|Iess| 



et al. 2012 ). Convection in the ice I shell is expected for 
small grain sizes ( |Barr and McKinnon , 2QQ7| [Mitri and| 



Showman 2008), but the long- wavelength topography 
of Titan implies shell thickness variations that convec- 
tion might eliminate, indicating that Titan's ice I shell 



may be currently conductive (Nimmo and Bills, 2010). 
Possibly this shell is at least in part clathrate hydrate, 
which has different transport properties. If the man- 
tle of high-pressure ice is convecting, then a system of 
equations, for the convective instability of and temper- 
ature contrasts across the thermal boundary layers at 
the interface with the ice/rock interior, is necessary to 
precisely calculate the thermal evolution of Titan. Be- 
cause our focus is on the thermal evolution of Titan's 
deep interior, we assume for simplicity that T5 is con- 
stant in all simulations. 



3. Numerical models 

Simulations are performed to calculate the thermal 
history of Titan's ice/rock interior for 4.5 Gyr. The 
model described above was numerically iterated with a 
time step of 1 Myr. Internal temperatures are compared 
to the phase diagram for pure water ice from |Petrenko| 
and Whitworth (1999) at each time step to determine 
whether melting occurs. Titan, as a whole, has p = 
1,881 kg m -3 and S « 0.44. We assume that Titan is 
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initially composed of an outer shell with a thickness of 
575 km and a density of 1,200 kg m -3 . Therefore, Ri = 
2,000 km and S ~ 0.82 in the ice/rock interior. The av- 
erage density of the ice/rock interior, p = 2,653 kg m -3 . 
In this simple two-layer model, the normalized moment 
of inertia coefficient C « 0.343. The melting tempera- 
tures of water ice are ~280 and 500 K at pressures of 
~0.7 and 4.7 GPa, respectively, which are the condi- 
tions at the top and bottom of the ice/rock interior. 

A wide parameter space was explored by varying the 
initial conditions for hundreds of simulations. The tem- 
perature at the top of the thermal boundary layer in the 
ice/rock interior was fixed at X5 = 200 K, but the effects 
of different initial temperature gradients were explored 
by using both T m (0) = 200 and 125 K, as suggested 



by models of Titan's accretion (Mitri et al. , 2010). Ini- 



tial rock mass fraction gradients were imposed that pro- 
duced effective density changes throughout the ice/rock 
interior as large as Ap s ~ O.lp (i.e., ~10%). An initial 
temperature gradient, if present, also produces an ef- 
fective density change of Ap t ~ aAT^, where AT^ = 
T5 — T m (0). In some simulations, no initial composi- 
tional gradient was imposed. Because of the large un- 
certainty on the viscosity of ice/rock mixtures, we tested 
effective viscosities near melting in the range v$ = 10 11 
to 10 15 Pa s. (Even larger values increase the likeli- 
hood of melting and lower values are much lower than 
physically plausible.) To account for the possibility of 
potassium leaching and uncertainties in the rock min- 
eralogy, we conducted simulations using three different 
depletion factors: / = 1.0, 0.7, and 0.3. Guided by the 
results of these simulations, we construct some plausible 
present-day internal structure models for Titan. 



4. Results 

The following sections summarize the results of our 
simulations, beginning with the presentation of four 
representative thermal histories. Then, we discuss the 
effects of various initial conditions using a large num- 
ber of simulation results. Finally, simulations in which 
melting occurred are analyzed to determine the details 
of melting in the ice/rock interior. 



Comparisons between pairs of these simulations illu- 
minate the effects of varying the magnitude of compo- 
sitional gradients, potassium depletion, and initial tem- 
perature gradients. Simulations shown in panels A and 
B, for instance, both have / = 1.0, v$ = 10 14 Pa s, and 
Tm(0) = 125 K. Panel A features a ~1% compositional 
gradient, i.e., Ap s = O.Olp, whereas the initial density 
gradient in panel B is only the result of the initial tem- 
perature gradient. The growth of the convecting layer 
is relatively delayed by ~500 Myr in panel A, allowing 
an additional ~100 K of heating in the ice/rock interior. 

Panels A and C illustrate the effects of potassium de- 
pletion on the thermal evolution of an ice/rock interior 
with a ~1% compositional gradient. All initial condi- 
tions are identical for both simulations, except / = 0.3 
in panel C. The primary effect of potassium depletion is 
to slow the growth of the convecting layer because heat- 
ing in the underlying, non-convecting layer is necessary 
to overcome the initial compositional gradient. Specifi- 
cally, convection in the entire ice/rock interior is delayed 
by ~1.25 Gyr. The maximum temperatures reached in 
the non-convecting ice/rock interior, however, are only 
slightly lower in panel C than in panel A. On the other 
hand, the temperature of the convecting layer is gen- 
erally reduced in panel C, since the convecting layer 
is relatively thin; internal heating is lessened; and the 
incorporation of hot material from the non-convecting 
interior is more gradual. 

Finally, panels C and D reveal the effects of varying 
the effective viscosity near melting. Both simulations 
include substantial potassium depletion (/ = 0.3) and 
start with a ~1% compositional gradient. While panel 
C features v$ = 10 14 Pa s as usual, vq = 10 12 Pa s in 
panel D. With a lower viscosity, the convecting layer 
has a lower temperature and takes slightly longer to 
encompass the entire interior. In fact, convecting layer 
growth is faster for the first ~500 Myr in panel D, but 
marginally slower thereafter. Thus, heat flux out of 
the convecting layer is relatively low after a period of 
relatively rapidly cooling and convecting layer growth. 
Compared to the other initial conditions, however, the 
sensitivity of the model to uncertainties in viscosity is 
very small. 



4.1. Sample thermal histories 

Four examples of calculated thermal histories are 
shown in Fig. [3j In each simulation, a convecting layer 
with constant temperature grows over time, eventually 
encompassing the entire ice/rock interior. From 2.5 
to 4.5 Gyr, this convecting ice/rock layer cools mono- 
tonically. Underneath the convecting layer, tempera- 
tures increase from radiogenic heating, although the 
initial linear temperature gradient is preserved. No 
melting occurs in any of these four examples within 
4.5 Gyr because the convecting layer grows fast enough 
to encompass the entire interior before temperatures in 
the non-convecting layer increase sufficiently to cause 
widespread melting and differentiation. 



4.2. Sensitivity to initial conditions 

Figure [4] shows results from 810 simulations of the 
thermal evolution of Titan's ice/rock interior. The pur- 
pose of each panel is to illustrate whether or not melting 
occurs for a particular set of initial conditions. Horizon- 
tal axes indicate the range of density differences across 
the ice/rock interior induced by an initial compositional 
gradient, i.e., Ap s = to 0.03p. Each vertical axis rep- 
resents a range of effective viscosities near melting, i.e., 
vo = 10 11 to 10 15 Pa s. Blue squares signify simulations 
in which melting never occurred, whereas red squares 
represent simulations where widespread melting would 
have caused differentiation of Titan's ice/rock interior 
within 4.5 Gyr. 
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Time [Myr] Time [Myr] 




Time [Myr] Time [Myr] 

Figure 3: Representative results for simulations in which melting did not occur. Temperature profiles in Titan's ice/rock interior are 
plotted until the time after which the convecting interior simply cools monotonically. Default initial conditions for the four simulations 
are T b = 200 K, T m (0) = 125 K, / = 1, log 10 (i/ ) = 14, and Ap s = O.Olp. Panel B has Ap s = 0. Panel C has / = 0.3. Panel D has / 
= 0.3 and log 10 (z/ ) = 12. 



Each panel represents a series of simulations that 
started with different initial temperature gradients and 
magnitudes of potassium depletion. The top and bot- 
tom rows of panels, for instance, have T m (0) = 125 and 
200 K, respectively. In general, melting is less likely 
to occur without an initial temperature gradient, be- 
cause a temperature decrease with depth induces a den- 
sity gradient that suppresses convection. High- viscosity 
simulations with low compositional gradients in pan- 
els D and E feature melting, while the corresponding 
simulations in panels A and B do not. The lack of a 
significant density gradient from compositional or tem- 
perature changes in these simulations causes the con- 
vecting layer to very quickly encompass the interior, 
homogenizing internal temperatures. Because the aver- 
age internal temperature is comparatively high at first, 
however, melting occurs near the top of the ice/rock 
interior, where melting temperatures are relatively low. 

Figure [4] also elucidates the effects of potassium de- 
pletion. From left to right, columns of panels have / = 
1.0 (A and D), 0.7 (B and E), and 0.3 (C and F). Reduc- 
ing the abundance of potassium in the rock component 
of the ice/rock interior only makes melting marginally 
less likely. As seen in the example thermal histories, 
potassium depletion only lengthens the time scales on 
which temperatures increase; the maximum tempera- 



tures reached in ice/rock interior are barely lessened. 
For any plausible choice of initial conditions, melting 
occurs if the initial compositional gradient is larger than 
-2%. 

4-3. When and where melting occurs 

Figure [5] contains more information about 173 sim- 
ulations in which melting occurred for initial compo- 
sitional gradients ranging from 1 to 10%. Specifically, 
panel A illustrates the time at which melting occurred 
and panel B shows the value of R C1 the radius of the non- 
convecting ice/rock interior, at the time when melting 
began. As in Fig.|4j six different sets of initial conditions 
were used. That is, black and red symbols respectively 
represent T m (0) = 125 and 200 K. For circular, square, 
and triangular markers, respectively, / = 1.0, 0.7, and 
0.3. With large density gradients, the non-convecting 
ice/rock core always melts before the ice in the convect- 
ing layer. 

The magnitudes of potassium depletion and the ini- 
tial temperature gradient more strongly affect when 
melting occurs than where it begins. That is, increasing 
potassium depletion can delay melting by ~500 Myr to 
1.75 Gyr for a given initial compositional gradient, but 
the radius of the non-convecting core at the time when 
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Figure 4: Sensitivity analyses results for 810 simulations exploring the effects of potassium depletion, effective viscosity near melting, and 
initial temperature and compositional gradients. In each panel, the horizontal axis shows the total density variation across the ice/rock 
interior produced by the initial compositional gradient and the vertical axis shows the dynamic viscosity near melting assumed for each 
simulation. Blue circles represent simulations in which no melting occurred in 4.5 Gyr; red squares signify simulations in which melting 
occurred at any point in the ice/rock interior. All simulations have Tb = 200 K. From left to right, each column of panels features / 
= 1.0, 0.7, and 0.3. The top and bottom rows of panels have T m (0) = 125 and 200 K, respectively. For initial compositional gradients 
producing Ap s > 0.02p, melting occurs during all simulations. 



melting begins only slightly decreases. With no ini- 
tial temperature gradient but a higher average internal 
temperature, melting occurs sooner, but the convecting 
layer has grown further before melting begins. As com- 
positional gradients increase, melting occurs sooner, the 
convecting layer is thinner when melting occurs, and 
simulation sensitivity to the choice of initial conditions 
decreases. 



5. Discussion 

5.1. Stability of a partially differentiated Titan 

A partially differentiated Titan is stable over geologic 
time if it can safeguard the ice in its deep interior from 
widespread melting caused by radiogenic heating. If 
thermal convection is vigorous throughout the entire 
ice/rock layer, then melting can be avoided. Even for 
small compositional gradients (~1-1.5%) a convecting 
layer can grow to encompass the entire ice/rock inte- 
rior before internal temperature profiles cross the melt- 
ing curve of water ice. Larger compositional gradients 
(~5-10%), however, are plausible since Titan likely ac- 
creted from heterogenous planetesimals. Thermal buoy- 
ancy simply cannot overcome these stabilizing gradi- 
ents without causing ice melting and complete differ- 
entiation. The likely presence of ammonia in the ice 
component, which may significantly lower the melting 



temperature (e.g., Loveday et al. 2009), compounds the 
problem. 

No plausible selection of initial conditions allows suf- 
ficiently efficient heat removal in the presence of a large 
compositional gradient. A substantial amount of heat- 
producing 40 K may be leached out of the rock com- 
ponent through contact with liquid water or ammonia- 
water solutions. Roughly 30% of the original abundance 
of potassium has been estimated to have been extracted 
from the silicates of Ganyme de (|Kirk and Stevenson| 



1987) and Titan (Engel et al. , 1994), and almost com- 



plete leaching of several elements has been suggested 
for Enceladus (iGlein and Shock 2010). A reduction 



in the magnitude of radiogenic heating may be impor- 
tant to preventing Titan's interior from reaching tem- 



peratures associated with silicate dehydration ( Cas tillo- 
Rogez and Lunine 2010| ), but, as seen in Fig. [2j any 
plausible abundance of radioisotopes, even with ^70% 
leaching of 40 K, will eventually produce enough heat 
to melt ice. Internal heating drives the growth of the 
double-diffusive convecting layer, so potassium deple- 
tion suppresses the initialization of convection through- 
out the entire ice/rock interior. 

Double-diffusive convection involves opposing gradi- 
ents of two components that have different diffusivities. 
With components like salt and water, both diffusivities 
are non-zero, and multiple convecting layers form. In 
Titan's ice/rock system, however, the diffusivity of the 
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Figure 5: Details of melting during 186 simulations of Titan's 
thermal evolution with initial compositional gradients that pro- 
duce density contrasts from O.Olp to O.lp across the ice/rock 
interior. Panel A shows when melting occurs in the non- 
convecting ice/rock interior. Panel B illustrates the size of the 
non-convecting ice/rock core at the time when melting occurs. 
All simulations have T& = 200 K. Black and red symbols respec- 
tively represent T m (0) = 125 and 200 K. Circular, square, and 
triangular markers represent / = 1.0, 0.7, and 0.3, respectively. 



rock component is effectively zero, and we only model 
the growth of a single layer. For small compositional 
gradients, the time scale for the growth of a convec- 
tive instability of a new layer is much greater the time 
scale on which the first layer grows. With sufficiently 
large compositional gradients, growth of the convecting 
layer is slow enough that additional layers could form if 
melting never occurred, but melting typically occurs in 
~250 to 750 Myr. In any case, heat transport through 
a series of convecting layers is comparatively inefficient 
(e.| 



Turner, 1974), and the formation of many layers 



is unlikely to prevent complete differentiation. 



5.2. Alternative internal structures 

A partially differentiated Titan is a non-unique in- 



terpretation of the available gravity data (less et al. 
2010). Because Titan's ice/rock interior is unstable 



to differentiation over geologic time, alternate models 
should be considered. One popular hypothesis is that 
Titan's interior is chiefly composed of hydrated silicates 
(e-g-, 



2010 



Fortes et al. 
Fortesl 120121). 



2007 Castillo-Rogez and Lunine 



Three serpentine polymorphs are 



lizardite, and antigorite. Antigorite is the most abun- 
dant, and its properties are usually used to track the 
thermochemical evolution of model assemblages (e.g., 
Castillo-Rogez and Lunine] |2010[ ). After heating melts 



the ice component of an ice/rock mixture, redistribut- 
ing rock and melt into a stable density structure with 
homogenized internal temperatures takes a few hundred 



Myr ( |Kirk and Stevenson| |1987 Lunine and Stevenson 



1987). The serpentinization reaction between the lib- 



erated silicates and the surrounding liquid would likely 
proceed to completion well before core overturn is fin- 



ished (MacDonald and Fyfe 1985). 



Models assuming a core composed entirely of antig- 
orite are missing a significant mass of iron. The for- 
mula for the magnesium end-member of antigorite is 
Mg3Si 2 5 (OH)4. But Fe 2+ and Fe 3+ often substitute 
into the Mg 2+ and Al 3+ sites, respectively, yielding a 
typic al iron number of Fe/(Fe+Mg) ~ 0.2 ( Scott et al. , 
|2002[ ). That is, antigorite has <14 wt% Fe, significantly 
less than the ~18.2 wt% Fe typically found in CI chon- 



drites (Hutchinson, 2004, p. 29). So, -4 wt% Fe, equiv- 



alent to an Fe-rich core with a radius of a few hundred 
km, must accompany a mass of antigorite created from 
the serpentinization of silicates with chondritic compo- 
sition. Some iron may be in soluble form in Titan's 
ocean, and the presence of a relatively small Fe-rich core 
may not dramatically affect Titan's moment of inertia 



(Castillo-Rogez and Lunine 2010]). 



Figure [6] shows two end-member models for the 
present-day internal structure of Titan for which C ~ 
0.34. In structure A, Titan is dominated by an hydrated 
silicate core with R = 1,950 km and p — 2,700 kg m -3 . 
An Fe-rich core with R = 500 km and p = 6,500 kg m -3 , 
representing ~4 wt% of the iron/rock interior is also 
present, although the formation of an Fe-rich core would 
likely dehydrate the overlying silicates. The overlying 
shell, composed of an ice I lid, an ocean, and high- 
pressure ice polymorphs, is modeled as a single layer 
with an average density of 1,202 kg m -3 . The lower 
~300 km of antigorite could be replaced with dehy- 
drated silicates in structure A without decreasing the 
Mol coefficient below C ~ 0.332 

A key question is whether significant dehydration of 
the silicate core should be expected. The equilibrium 
temperature for the antigorite dehydration reaction is 
^800 to 900 K for conditions appropriate to Titan (e.g., 



Perrillat et al. , 2005), although dehydration in a water 



saturated environment may be kinetically hindered and 
not occur until temperatures reach ~1060 K (Seipold 



and Schilling 2003 ) , and the heat flux from the advec- 



tion of warm water could be significant. Castillo-R ogez 
and Lunine ( |2010[ ) found that conductive heat trans- 



port could delay significant dehydration until very re- 
cently, assuming ~30% potassium depletion and a rela- 
tively high specific heat for antigorite. Serpentine dehy- 
dration at the present may be responsible for melting 
clathrate hydrates in the ov erlying ice shell, releasing 
methane (Tobie et al. 2006) and radiogenic 40 Ar (Nie- 



typically found in carbonaceous chondrites: chrysotile, |mann et al. 2005) onto the surface and into the at 
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Figure 6: Two end-member internal structures for Titan that have C ~ 0.34. Possible material layers include an Fe-rich core, dehydrated 
silicates, hydrated silicates, and an ice/rock mixture. Each model includes high-pressure ices, a water ocean with dissolved gases and 
solids, and a layer of water ice I, all of which are modeled as a single shell with constant density. Titan's atmosphere is neither depicted 
nor modeled. 



mosphere. We suggest that, if Titan formed partially 
differentiated, double-diffusive convection in the initial 
ice/rock mixture may have delayed core overturn for 
an additional ~1 Gyr beyond what was assumed in 



Castillo-Rogez and Lunine (2010), helping the eventual 



core of hydrated silicates avoid dehydration. 

Current data do not exclude partially differentiated 
models featuring a large silicate core. Structure B in 
Fig. [6] is an example of such a model, similar to those 
considered in McKinnon and Bland (2011 ). An ice rock 
mantle with R = 2,075 km, p = 2,200 kg m~ 3 , and S 
~ 0.63 overlies a core of dehydrated silicates with R 
= 1,300 km and p = 3,300 kg m -3 . A decreased rock 
mass fraction and a smaller length scale for the ice/rock 
layer might allow efficient convection, even in the pres- 
ence of a stabilizing compositional gradient. Evaluat- 
ing the thermal stability of such a structure would re- 
quire coupling the thermal evolution of all three lay- 
ers, potentially including the effects of tidal dissipa- 



tion, as has been done for Ganymede (Bland et al. 



2009). More complicated variations on structure B fea- 



turing different ice/rock mixtures or a hydrated por- 
tion of the silicate core are also valid interpretations of 
the gravity data. Differences in the predicted tidal re- 
sponses of these two models may allow discrimination 
between these two alternatives ( Castillo-Rogez and Lu 



2012]). 



5.3. Implications for other icy satellites 

Gravity measurements by spacecraft suggest a spec- 
trum of internal structures for the icy satellites of 
Jupiter and Saturn. The formation of an undifferenti- 
ated or a partially differentiated satellite requires a se- 
quence of assumptions: slow accretion in a gas-starved 



disk, delayed to avoid catastrophic heating from the ra- 
dioactive decay of short-lived isotopes, and then sur- 
vival of intense impacts from the Late Heavy Bombard- 
ment ( |Barr and Canup[ |2008| |2010| |Barr et ah} |2010[ ). 
The Radau-Darwin approximation is primarily used to 
calculate the Mol of a satellite with the measured J2 
and C22 coefficients, but a unique solution only exists 
for a satellite in hydrostatic equilibrium. 



Non-hydrostatic structures in the degree 2 gravity 
of Titan may introduce significant inaccuracies into 



the Radau-Darwin approximatio n (|Gao and Steven 



son 



2012). Indeed, less et al. (2010) reported non- 



hydrostatic geoid height variations as large as 19 m 
for Titan, possibly indicating that Titan's Mol coef- 
ficient could closer to that of Ganymede than earlier 
reported. The measured Mol for Callisto is likewise 
subject to large uncertainties because of possible non- 
hydrostatic effects, which are not excluded by current 



data (Anderson et al. , 2001), but an equivalent depar- 



ture from hydrostatic equilibrium would have a rela- 
tively small effect on Ganymede. The Galilean satel- 
lites are not as compositionally heterogeneous as the 
major Saturnian satellites (Anderson et al. 1996b| 



1998| |2001[ ), so the likelihood of substantial composi- 



tional gradients suppressing convection in Callisto may 
be lessened. In any case, the radio science instrument 
on the planned JUICE (JUpiter Icy moons Explorer) 
mission may allow the vastly improved measurement of 
the gravity fields of Ganymede and Callisto and thus 
the identification of gravity anomalies and significant 
non-hydrostatic features. With these new data, we may 
continue to vet the spectrum of internal differentiation 
states proposed for icy satellites. 
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6. Conclusions 

Titan can survive accretion and the Late Heavy Bom- 
bardment without suffering complete differentiation if it 
formed slowly in a gas-starved disk, and a partially dif- 
ferentiated Titan is consistent with recently obtained 
gravity data. But radiogenic heating tends to melt 
ice in the undifferentiated interior over geologic time, 
causing irreversible sinking of the rock component. If 
the ice/rock interior is initially homogenous, then ther- 
mal convection may remove sufficient heat to prevent 
melting. Two sources of compositional gradients, how- 
ever, act to suppress convection. First, Titan likely 
had a temperature profile immediately after accretion in 
which temperature decreased with depth. Second, ac- 
cretion from heterogenous planetesimals would produce 
a stabilizing rock mass fraction gradient, where rock 
mass fraction and thus density increased with depth. 
Convection in the presence of opposing compositional 
and thermal gradients is "double-diffusive", compris- 
ing a convecting layer that slowly grows from the top 
of the ice/rock interior. For compositional gradients 
larger than a few percent, this convecting layer does 
not grow fast enough to prevent widespread melting in 
the non-convecting core. A partially differentiated Ti- 
tan is therefore unstable over geologic time. The general 
instability of ice/rock mixtures with compositional gra- 
dients may have significant consequences for other icy 
satellites. 
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